;********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;load library for  wetbulb_stull
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/heat_stress.ncl"
;The *area_conserve_remap *and *area_conserve_remap_Wrap* have been *deprecated*.
;The* ESMF regrid* <http://www.ncl.ucar.edu/Applications/ESMF.shtml> library is the recommended replacement. It uses a superior methodology.
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
;********************************************
 begin

;set baseline
 dirPath = "~/source_data/"
 dirPath2 = dirPath+"/outdoor/"

 mo  = (/"ACCESS-CM2","AWI-CM-1-1-MR","BCC-CSM2-MR","CMCC-CM2-SR5","EC-Earth3","GFDL-ESM4","IITM-ESM", \
         "KIOST-ESM","MIROC6","MIROC-ES2L","MPI-ESM1-2-HR","MRI-ESM2-0"/); 12 models
 case = (/"ssp585" /)
 ncase = dimsizes(case);
 nmod = dimsizes(mo)

 syr = 1990; including historical data to identify 1C warming
 fyr = 2099;
 csyr = sprinti("%0.4i",syr);
 cfyr = sprinti("%0.4i",fyr);
 nyr = fyr - syr + 1;
 ndays = nyr*365

;read dimension for 1 deg data
  dimPath=dirPath+"/domain/"
  landfrac = "sftlf_fx_360x180_1deg_gn.nc"
  area = "areacella_fx_360x180_1deg_gn.nc"

  fdim  = addfile(dimPath+landfrac, "r")
  fdim2 = addfile(dimPath+area, "r")
  lat := fdim->lat
  lon := fdim->lon
  nlat = dimsizes(lat);
  nlon = dimsizes(lon);
  garea = fdim2->areacella  ;m^2
  lfrac = fdim->sftlf/100; %
  wgt := garea*1e-6*lfrac ; sqrkm
  wgt@_FillValue = lfrac@_FillValue; make sure FillValue is set correctly for wgt
  copy_VarCoords(garea, wgt)
  printVarSummary(wgt)
  print("global total land area: "+sum(wgt)+" sq km")



  ;----data from Data.Outdoor_Gday_impacts_gridcell-ens-stat-vcbc-MH300.ncl
  ;spatial statistics and extent of impact
  ;*masked by water access data
  csv_file1 = dirPath2+"Gday0_area_pop_sun_warminglevels_10yravg_ens_median-vcbc-outdoor-water.csv"
  csv_file2 = dirPath2+"Gday0_area_pop_diff_warminglevels_10yravg_ens_median-vcbc-outdoor-water.csv"
  csv_file3 = dirPath2+"Gday0_area_pop_zero_warminglevels_10yravg_ens_median-vcbc-outdoor-water.csv"

  ;masked by rural population data
  csv_file4 = dirPath2+"Gday0_area_pop_sun_warminglevels_10yravg_ens_median-vcbc-outdoor-rural.csv"
  csv_file5 = dirPath2+"Gday0_area_pop_diff_warminglevels_10yravg_ens_median-vcbc-outdoor-rural.csv"
  csv_file6 = dirPath2+"Gday0_area_pop_zero_warminglevels_10yravg_ens_median-vcbc-outdoor-rural.csv"


;---Read in file as array of strings so we can parse each line
   lines1 := asciiread(csv_file1,-1,"string")
   lines2 := asciiread(csv_file2,-1,"string")
   lines3 := asciiread(csv_file3,-1,"string")
   lines4 := asciiread(csv_file4,-1,"string")
   lines5 := asciiread(csv_file5,-1,"string")
   lines6 := asciiread(csv_file6,-1,"string")
   nlines = dimsizes(lines1)-1   ; First line is a header
   delim = ","

  ;---skip the header line
   lines1 := lines1(1:)
   lines2 := lines2(1:)
   lines3 := lines3(1:)
   lines4 := lines4(1:)
   lines5 := lines5(1:)
   lines6 := lines6(1:)

;---Read fields (columns)
   levs =  tofloat(str_get_field(lines1,1,delim))
   ;--select the warming levels to plot
   ;levels := fspan(1,4,7); 1 to 4 degC by 0.5 degC
   levels := fspan(1,4,31); 1 to 4 degC by 0.1 degC
   levid = get1Dindex(levs,levels)
   nlev = dimsizes(levels)

   ;area and population impact (*masked by water access)
   G0_warea_sun_median =  tofloat(str_get_field(lines1(levid),2,delim))
   G0_wpop_sun_median  =  tofloat(str_get_field(lines1(levid),3,delim))
   G0_warea_sun_high =  tofloat(str_get_field(lines1(levid),4,delim))
   G0_wpop_sun_high  =  tofloat(str_get_field(lines1(levid),5,delim))
   G0_warea_sun_low =  tofloat(str_get_field(lines1(levid),6,delim))
   G0_wpop_sun_low  =  tofloat(str_get_field(lines1(levid),7,delim))

   G0_warea_diff_median =  tofloat(str_get_field(lines2(levid),2,delim))
   G0_wpop_diff_median  =  tofloat(str_get_field(lines2(levid),3,delim))
   G0_warea_diff_high =  tofloat(str_get_field(lines2(levid),4,delim))
   G0_wpop_diff_high  =  tofloat(str_get_field(lines2(levid),5,delim))
   G0_warea_diff_low =  tofloat(str_get_field(lines2(levid),6,delim))
   G0_wpop_diff_low  =  tofloat(str_get_field(lines2(levid),7,delim))
   G0_warea_diff1_median =  tofloat(str_get_field(lines2(levid),8,delim))
   G0_wpop_diff1_median  =  tofloat(str_get_field(lines2(levid),9,delim))
   G0_warea_diff2_median =  tofloat(str_get_field(lines2(levid),10,delim))
   G0_wpop_diff2_median  =  tofloat(str_get_field(lines2(levid),11,delim))

   G0_warea_zero_median =  tofloat(str_get_field(lines3(levid),2,delim))
   G0_wpop_zero_median  =  tofloat(str_get_field(lines3(levid),3,delim))
   G0_warea_zero_high =  tofloat(str_get_field(lines3(levid),4,delim))
   G0_wpop_zero_high  =  tofloat(str_get_field(lines3(levid),5,delim))
   G0_warea_zero_low =  tofloat(str_get_field(lines3(levid),6,delim))
   G0_wpop_zero_low  =  tofloat(str_get_field(lines3(levid),7,delim))

   ;area and population impact (*masked by rural population data)
   G0_rarea_sun_median =  tofloat(str_get_field(lines4(levid),2,delim))
   G0_rpop_sun_median  =  tofloat(str_get_field(lines4(levid),3,delim))
   G0_rarea_sun_high =  tofloat(str_get_field(lines4(levid),4,delim))
   G0_rpop_sun_high  =  tofloat(str_get_field(lines4(levid),5,delim))
   G0_rarea_sun_low =  tofloat(str_get_field(lines4(levid),6,delim))
   G0_rpop_sun_low  =  tofloat(str_get_field(lines4(levid),7,delim))

   G0_rarea_diff_median =  tofloat(str_get_field(lines5(levid),2,delim))
   G0_rpop_diff_median  =  tofloat(str_get_field(lines5(levid),3,delim))
   G0_rarea_diff_high =  tofloat(str_get_field(lines5(levid),4,delim))
   G0_rpop_diff_high  =  tofloat(str_get_field(lines5(levid),5,delim))
   G0_rarea_diff_low =  tofloat(str_get_field(lines5(levid),6,delim))
   G0_rpop_diff_low  =  tofloat(str_get_field(lines5(levid),7,delim))
   G0_rarea_diff1_median =  tofloat(str_get_field(lines5(levid),8,delim))
   G0_rpop_diff1_median  =  tofloat(str_get_field(lines5(levid),9,delim))
   G0_rarea_diff2_median =  tofloat(str_get_field(lines5(levid),10,delim))
   G0_rpop_diff2_median  =  tofloat(str_get_field(lines5(levid),11,delim))

   G0_rarea_zero_median =  tofloat(str_get_field(lines6(levid),2,delim))
   G0_rpop_zero_median  =  tofloat(str_get_field(lines6(levid),3,delim))
   G0_rarea_zero_high =  tofloat(str_get_field(lines6(levid),4,delim))
   G0_rpop_zero_high  =  tofloat(str_get_field(lines6(levid),5,delim))
   G0_rarea_zero_low =  tofloat(str_get_field(lines6(levid),6,delim))
   G0_rpop_zero_low  =  tofloat(str_get_field(lines6(levid),7,delim))



  filNc1  = "water_pop_warminglevels_12models-ens-median.nc"
  pthNc1  = dirPath2+filNc1
  ncdf1 = addfile(pthNc1 ,"r")        ; open output netCDF file
  WP_lev_med = ncdf1->WP_lev_med
  level1 = ncdf1->lev; 0.8 to 4.5 by 0.1 degC
  WP_lev = dim_sum_n(WP_lev_med, (/1,2/))*1e-6
  print("")
  print("")
  print("population without access to piped water at warming level:"+level1(2)+" : "+WP_lev(2)+" million")

  filNc2  = "rural_pop_warminglevels_12models-ens-median.nc"
  pthNc2  = dirPath2+filNc2
  ncdf2 = addfile(pthNc2 ,"r")        ; open output netCDF file
  RP_lev_med = ncdf2->RP_lev_med
  level2 = ncdf2->lev; 0.8 to 4.5 by 0.1 degC
  RP_lev = dim_sum_n(RP_lev_med, (/1,2/))*1e-6
  print("rural population at warming level:"+level2(2)+" : "+RP_lev(2)+" million")
  levid2 = get1Dindex(level2, levels); select levels 1-4 by 0.1 for the plot

  
  levs1 = ispan(1,4,1)
  levid1 = get1Dindex(levels, levs1); select levels 1-4 by 1 for summary statistics
  print("")
  print("summary data used in the main text")
  print("")

  ;mask global land area by water access data and rural population data
  water_area = where(WP_lev_med(0,:,:) .gt. 0, wgt, wgt@_FillValue); water area 
  rural_area = where(RP_lev_med(0,:,:) .gt. 0, wgt, wgt@_FillValue); rural area 
  TOTWA = sum(water_area)*100*1e-6; Mha
  TOTRA = sum(rural_area)*100*1e-6; Mha
  print(" ")
  print("global total area without indoor water access in 2022: "+TOTWA+ " Mha")
  print("global total rural area in 2020: "+TOTRA+ " Mha")

  print(" ")

  ;fraction of impacted area
  G0_wareafrc_zero_median := G0_warea_zero_median/TOTWA*100
  G0_wareafrc_diff_median := G0_warea_diff_median/TOTWA*100
  G0_wareafrc_diff1_median := G0_warea_diff1_median/TOTWA*100
  G0_wareafrc_diff2_median := G0_warea_diff2_median/TOTWA*100
  G0_wareafrc_sun_median := G0_warea_sun_median/TOTWA*100
  G0_wareafrc_zero_high := G0_warea_zero_high/TOTWA*100
  G0_wareafrc_diff_high := G0_warea_diff_high/TOTWA*100
  G0_wareafrc_sun_high := G0_warea_sun_high/TOTWA*100
  G0_wareafrc_zero_low := G0_warea_zero_low/TOTWA*100
  G0_wareafrc_diff_low := G0_warea_diff_low/TOTWA*100
  G0_wareafrc_sun_low := G0_warea_sun_low/TOTWA*100

  G0_rareafrc_zero_median := G0_rarea_zero_median/TOTRA*100
  G0_rareafrc_diff_median := G0_rarea_diff_median/TOTRA*100
  G0_rareafrc_diff1_median := G0_rarea_diff1_median/TOTRA*100
  G0_rareafrc_diff2_median := G0_rarea_diff2_median/TOTRA*100
  G0_rareafrc_sun_median := G0_rarea_sun_median/TOTRA*100
  G0_rareafrc_zero_high := G0_rarea_zero_high/TOTRA*100
  G0_rareafrc_diff_high := G0_rarea_diff_high/TOTRA*100
  G0_rareafrc_sun_high := G0_rarea_sun_high/TOTRA*100
  G0_rareafrc_zero_low := G0_rarea_zero_low/TOTRA*100
  G0_rareafrc_diff_low := G0_rarea_diff_low/TOTRA*100
  G0_rareafrc_sun_low := G0_rarea_sun_low/TOTRA*100


  print(" ")
  print("water area with Gday >0")
  print("under zero solar radiation:")
  print("per warming level:"+levels(levid1)+" : "+G0_warea_zero_median(levid1)+" Mha; "+G0_wareafrc_zero_median(levid1)+"%")
  print("under diffuse radiation:")
  print("per warming level:"+levels(levid1)+" : "+G0_warea_diff_median(levid1)+" Mha; "+G0_wareafrc_diff_median(levid1)+"%")
  print("-50% of diffuse, per warming level:"+levels(levid1)+" : "+G0_warea_diff1_median(levid1)+" Mha (diffuse lower)")
  print("+50% of diffuse, per warming level:"+levels(levid1)+" : "+G0_warea_diff2_median(levid1)+" Mha (diffuse upper)")
  print("under full solar:")
  print("per warming level:"+levels(levid1)+" : "+G0_warea_sun_median(levid1)+" Mha; "+G0_wareafrc_sun_median(levid1)+"%")
  print("high percentile per warming level:"+levels(levid1)+" : "+G0_warea_sun_high(levid1)+" Mha (full solar) "+G0_wareafrc_sun_high(levid1)+"%")

  print(" ")
  print("rural area with Gday >0")
  print("under zero solar radiation:")
  print("per warming level:"+levels(levid1)+" : "+G0_rarea_zero_median(levid1)+" Mha; "+G0_rareafrc_zero_median(levid1)+"%")
  print("under diffuse radiation:")
  print("per warming level:"+levels(levid1)+" : "+G0_rarea_diff_median(levid1)+" Mha; "+G0_rareafrc_diff_median(levid1)+"%")
  print("-50% of diffuse, per warming level:"+levels(levid1)+" : "+G0_rarea_diff1_median(levid1)+" Mha (diffuse lower)")
  print("+50% of diffuse, per warming level:"+levels(levid1)+" : "+G0_rarea_diff2_median(levid1)+" Mha (diffuse upper)")
  print("under full solar:")
  print("per warming level:"+levels(levid1)+" : "+G0_rarea_sun_median(levid1)+" Mha; "+G0_rareafrc_sun_median(levid1)+"%")
  print("high percentile per warming level:"+levels(levid1)+" : "+G0_rarea_sun_high(levid1)+" Mha (full solar) "+G0_rareafrc_sun_high(levid1)+"%")



  ;---present % impact together population count---

  G0_wpct_zero_median := G0_wpop_zero_median/WP_lev(levid2)*100
  G0_wpct_diff_median := G0_wpop_diff_median/WP_lev(levid2)*100
  G0_wpct_diff1_median := G0_wpop_diff1_median/WP_lev(levid2)*100
  G0_wpct_diff2_median := G0_wpop_diff2_median/WP_lev(levid2)*100
  G0_wpct_sun_median := G0_wpop_sun_median/WP_lev(levid2)*100
  G0_wpct_zero_high := G0_wpop_zero_high/WP_lev(levid2)*100
  G0_wpct_diff_high := G0_wpop_diff_high/WP_lev(levid2)*100
  G0_wpct_sun_high := G0_wpop_sun_high/WP_lev(levid2)*100
  G0_wpct_zero_low := G0_wpop_zero_low/WP_lev(levid2)*100
  G0_wpct_diff_low := G0_wpop_diff_low/WP_lev(levid2)*100
  G0_wpct_sun_low := G0_wpop_sun_low/WP_lev(levid2)*100


  G0_rpct_zero_median := G0_rpop_zero_median/RP_lev(levid2)*100
  G0_rpct_diff_median := G0_rpop_diff_median/RP_lev(levid2)*100
  G0_rpct_diff1_median := G0_rpop_diff1_median/RP_lev(levid2)*100
  G0_rpct_diff2_median := G0_rpop_diff2_median/RP_lev(levid2)*100
  G0_rpct_sun_median := G0_rpop_sun_median/RP_lev(levid2)*100
  G0_rpct_zero_high := G0_rpop_zero_high/RP_lev(levid2)*100
  G0_rpct_diff_high := G0_rpop_diff_high/RP_lev(levid2)*100
  G0_rpct_sun_high := G0_rpop_sun_high/RP_lev(levid2)*100
  G0_rpct_zero_low := G0_rpop_zero_low/RP_lev(levid2)*100
  G0_rpct_diff_low := G0_rpop_diff_low/RP_lev(levid2)*100
  G0_rpct_sun_low := G0_rpop_sun_low/RP_lev(levid2)*100

  print(" ")
  print("% water population with Gday >0")
  print("under zero solar radiation:")
  print("per warming level:"+levels(levid1)+" : "+G0_wpop_zero_median(levid1)+" million; "+G0_wpct_zero_median(levid1)+"%")
  print("under diffuse radiation:")
  print("per warming level:"+levels(levid1)+" : "+G0_wpop_diff_median(levid1)+" million; "+G0_wpct_diff_median(levid1)+"%")
  print("-50% of diffuse, per warming level:"+levels(levid1)+" : "+G0_wpop_diff1_median(levid1)+" million (diffuse lower)")
  print("+50% of diffuse, per warming level:"+levels(levid1)+" : "+G0_wpop_diff2_median(levid1)+" million (diffuse upper)")
  print("under full solar:")
  print("per warming level:"+levels(levid1)+" : "+G0_wpop_sun_median(levid1)+" million; "+G0_wpct_sun_median(levid1)+"%")
  print("high percentile per warming level:"+levels(levid1)+" : "+G0_wpop_sun_high(levid1)+" million (full solar) "+G0_wpct_sun_high(levid1)+"%")

  print(" ")
  print("% rural population with Gday >0")
  print("under zero solar radiation:")
  print("per warming level:"+levels(levid1)+" : "+G0_rpop_zero_median(levid1)+" million; "+G0_rpct_zero_median(levid1)+"%")
  print("under diffuse radiation:")
  print("per warming level:"+levels(levid1)+" : "+G0_rpop_diff_median(levid1)+" million; "+G0_rpct_diff_median(levid1)+"%")
  print("-50% of diffuse, per warming level:"+levels(levid1)+" : "+G0_rpop_diff1_median(levid1)+" million (diffuse lower)")
  print("+50% of diffuse, per warming level:"+levels(levid1)+" : "+G0_rpop_diff2_median(levid1)+" million (diffuse upper)")
  print("under full solar:")
  print("per warming level:"+levels(levid1)+" : "+G0_rpop_sun_median(levid1)+" million; "+G0_rpct_sun_median(levid1)+"%")
  print("high percentile per warming level:"+levels(levid1)+" : "+G0_rpop_sun_high(levid1)+" million (full solar) "+G0_rpct_sun_high(levid1)+"%")

exit

  ;--------plot a,b for water collection; c,d for farming
   xx = levels

   xp    = new( (/2*nlev/), float )
   yp1   = new( (/2*nlev/), float )
   yp2   = new( (/2*nlev/), float )
   yp3   = new( (/2*nlev/), float )
   yp4   = new( (/2*nlev/), float )
   yp5   = new( (/2*nlev/), float )
   yp6   = new( (/2*nlev/), float )
   yp11   = new( (/2*nlev/), float )
   yp22   = new( (/2*nlev/), float )
   yp33   = new( (/2*nlev/), float )
   yp44   = new( (/2*nlev/), float )
   yp55   = new( (/2*nlev/), float )
   yp66   = new( (/2*nlev/), float )

   xp(0:(nlev-1)) = xx
   xp(nlev:(2*nlev-1)) = xx(::-1) ; reverse

  ;ensemble 25-75th percentile CI for area and population impact
   yp1(0:(nlev-1)) = G0_warea_sun_high ; uppper bound
   yp1(nlev:(2*nlev-1)) = G0_warea_sun_low(::-1) ; lower bound
   yp2(0:(nlev-1)) = G0_warea_diff_high ; uppper bound
   yp2(nlev:(2*nlev-1)) = G0_warea_diff_low(::-1) ; lower bound
   yp3(0:(nlev-1)) = G0_warea_zero_high ; uppper bound
   yp3(nlev:(2*nlev-1)) = G0_warea_zero_low(::-1) ; lower bound

   yp11(0:(nlev-1)) = G0_wpop_sun_high ; uppper bound
   yp11(nlev:(2*nlev-1)) = G0_wpop_sun_low(::-1) ; lower bound
   yp22(0:(nlev-1)) = G0_wpop_diff_high ; uppper bound
   yp22(nlev:(2*nlev-1)) = G0_wpop_diff_low(::-1) ; lower bound
   yp33(0:(nlev-1)) = G0_wpop_zero_high ; uppper bound
   yp33(nlev:(2*nlev-1)) = G0_wpop_zero_low(::-1) ; lower bound

   yp4(0:(nlev-1)) = G0_rarea_sun_high ; uppper bound
   yp4(nlev:(2*nlev-1)) = G0_rarea_sun_low(::-1) ; lower bound
   yp5(0:(nlev-1)) = G0_rarea_diff_high ; uppper bound
   yp5(nlev:(2*nlev-1)) = G0_rarea_diff_low(::-1) ; lower bound
   yp6(0:(nlev-1)) = G0_rarea_zero_high ; uppper bound
   yp6(nlev:(2*nlev-1)) = G0_rarea_zero_low(::-1) ; lower bound

   yp44(0:(nlev-1)) = G0_rpop_sun_high ; uppper bound
   yp44(nlev:(2*nlev-1)) = G0_rpop_sun_low(::-1) ; lower bound
   yp55(0:(nlev-1)) = G0_rpop_diff_high ; uppper bound
   yp55(nlev:(2*nlev-1)) = G0_rpop_diff_low(::-1) ; lower bound
   yp66(0:(nlev-1)) = G0_rpop_zero_high ; uppper bound
   yp66(nlev:(2*nlev-1)) = G0_rpop_zero_low(::-1) ; lower bound



 ;---draw line plots
 imagePath = dirPath+"/../image/" 
 image1 = imagePath+"/Fig.4-newRH"
 wks1 = gsn_open_wks("pdf",image1)

 res = True
 res@gsnDraw              = False
 res@gsnFrame             = False

 ;res@vpXF                  = 0.1;0.15
 ;res@vpWidthF              = 0.65
 ;res@vpHeightF             = 0.50
 res@vpWidthF             = 0.32; 0.6 ;plot width
 res@vpHeightF            = 0.3; 0.3   ;plot length

 res@tiXAxisFontHeightF   = 0.014;0.022
 res@tiYAxisFontHeightF   = 0.014;0.022
 res@tmXBLabelFontHeightF = 0.012
 res@tmYLLabelFontHeightF = 0.012
 res@tmYRLabelFontHeightF = 0.012
 res@gsnStringFontHeightF = 0.014
 res@tmXTOn               = False
 ;res@tmXBMinorOn          = False
 ;res@tmBorderThicknessF    = 1.5

 res@pmLegendDisplayMode    = "Never";"Always"
 res@pmLegendSide           = "Top"; reference axis
 res@pmLegendParallelPosF   = .38; left-right adjustment
 res@pmLegendOrthogonalPosF = -.31;
 res@pmLegendWidthF         = 0.4;0.1
 res@pmLegendHeightF        = 0.04;0.4
 res@lgLabelFontHeightF     = .015
 res@lgPerimOn              = False
 res@lgOrientation = "horizontal"
 res@lgLabelAngleF = 300

 cmap = read_colormap_file("srip_reanalysis");("radar")
 cmap = cmap(::-1,:) ; reverse color map
 cmap := cmap((/1,2,3,4,5,6,7,8,10,11,12,13,15,16,18,0/),:); first black, last gray for ensemble median
 res@tiXAxisString   = ""; "Warming since pre-industrial (~S~o~N~C)" 
 res@trXMinF          = 1
 res@trXMaxF          = 4

 ;res@tfPolyDrawOrder    = "Predraw"               ; put line on top


 ;----Fig4a,b focus on water collection -----
  plot1 = new(3,graphic)
  plot2 = new(3,graphic)

  res@vpYF          = 0.90   ; location of plot1
  res@vpXF          = 0.10    ; location of plot1

  res@tmYLOn              = True
  res@tmYLLabelsOn        = True
  res@tmXBOn              = True
  res@tmXBLabelsOn        = True
  res@trYAxisType = "LinearAxis"
  res@gsnStringFontHeightF = 0.014;
  res@gsnLeftString       = "~F22~a";
  res@gsnCenterString     = "Focus on people collecting water"; 
  res@gsnLeftStringOrthogonalPosF = 0.03 ;adjust the space/distance of gsnLeftString/gsnCenterString
  res@gsnCenterStringOrthogonalPosF = 0.02 ;adjust the space/distance of gsnLeftString/gsnCenterString
  res@tiYAxisString    = "Land area with G~S~day~N~>0 (Mha)"; 
  res@tiYAxisJust      = "CenterCenter"
  res@tiYAxisOffsetXF  = -0.01; -0.01
  res@trYMinF             = 0; 
  res@trYMaxF             = 2466; 2160; Mha
  res@tmYROn              = False
  res@tmYLMode            = "Manual"
  res@tmYLTickStartF      = 0
  res@tmYLTickEndF        = 2400; 2100
  res@tmYLTickSpacingF    = 500
  res@tmYLMinorPerMajor   = 4
  res@tmYLAutoPrecision   = False
  res@tmYLPrecision       = 4

 ;resources for "right" Y axis
  resR = True
  resR@vpYF          = 0.50   ; location of plot1
  resR@vpXF          = 0.10    ; location of plot1
  resR@tiYAxisString = "Percent area (%)";+" [dashed]"
  resR@trYMinF         = 0
  resR@trYMaxF         = 40; 35
  resR@xyDashPatterns    = 0; 
  resR@xyLineThicknesses  = 1.5
  resR@xyLineOpacities    = 1
  resR@pmLegendDisplayMode    = "Never"

  resR@tmYRMode            = "Manual"
  resR@tmYRTickStartF      = 0
  resR@tmYRTickEndF        = 40; 35
  resR@tmYRTickSpacingF    = 5
  resR@tmYRMinorPerMajor   = 4
  resR@tmYRAutoPrecision   = False
  resR@tmYRPrecision       = 2

  res@xyLineColors     = "purple2";
  resR@xyLineColors    =  res@xyLineColors
  res@xyDashPatterns   = 0
  ;plot1(0) = gsn_csm_xy(wks1, xx, G0_warea_sun_median, res); full solar
  plot1(0) = gsn_csm_xy2(wks1, xx, G0_warea_sun_median, G0_wareafrc_sun_median, res, resR); 
 
  gsres                   = True                        ; poly res
  gsres@gsFillOpacityF    = 0.1
  gsres@gsFillColor       = "purple"                    ; color chosen
  dummy2 = gsn_add_polygon (wks1,plot1(0),xp,yp1,gsres);
  res@xyDashPatterns   = 1

  res@tiXAxisString    = " "
  res@tiYAxisString    = " "
  res@gsnLeftString    = " ";
  res@gsnCenterString    = " ";
  res@tmXBOn              = False
  res@tmYLOn              = False
  res@tmYLLabelsOn        = False
  res@tmXBLabelsOn        = False
  resR@tmYROn              = False
  resR@tmYRLabelsOn        = False
  resR@tiYAxisString      = ""

  res@xyLineColors     = "orange2";
  res@xyDashPatterns   = 0; 1
  resR@xyLineColors    =  res@xyLineColors

  ;plot1(1) = gsn_csm_xy(wks1, xx, G0_warea_diff_median, res); diffuse
  plot1(1) = gsn_csm_xy2(wks1, xx, G0_warea_diff_median, G0_wareafrc_diff_median, res, resR); 
  gsres@gsFillColor       = "orange"                    ; color chosen
  dummy1 = gsn_add_polygon (wks1,plot1(1),xp,yp2,gsres);
  res@xyDashPatterns   = 1

  res@xyLineColors      = "slategray"
  resR@xyLineColors    =  res@xyLineColors
  res@xyDashPatterns    = 0; 3
  ;plot1(2) = gsn_csm_xy(wks1, xx, G0_warea_zero_median, res);
  plot1(2) = gsn_csm_xy2(wks1, xx, G0_warea_zero_median, G0_wareafrc_zero_median, res, resR); 
  gsres@gsFillColor       = "slategray4"               ; color chosen
  dummy3 = gsn_add_polygon (wks1,plot1(2),xp,yp3,gsres);

  print("finish plot a")


 ;use vpXF to place two plots horizontally
  res@vpXF          = 0.6       ; location of plot2
  resR@vpXF          = 0.6       ; location of plot2
  res@tmYLOn              = True
  res@tmYLLabelsOn        = True
  res@tmXBOn              = True
  res@tmXBLabelsOn        = True

  res@gsnLeftString       = "~F22~b";
  res@gsnCenterString     = "Focus on people collecting water"; "Global";
  res@gsnLeftStringOrthogonalPosF = 0.04 ;adjust the space/distance of gsnLeftString/gsnCenterString
  res@gsnCenterStringOrthogonalPosF = 0.03 ;adjust the space/distance of gsnLeftString/gsnCenterString
  ;res@tiXAxisString    = "Warming since pre-industrial (~S~o~N~C)"
  res@tiYAxisString    = "Population with G~S~day~N~>0 (million)"; 
  res@tiYAxisJust      = "CenterCenter"
  res@tiYAxisOffsetXF  = 0;0.01
  res@trYMinF             = 0; 
  res@trYMaxF             = 281; 250; million
  res@tmYLOn              = True
  res@tmYROn              = False
  res@tmYLLabelsOn        = True
  res@tmYRLabelsOn        = False
  res@tmYLMode            = "Manual"
  res@tmYLTickStartF      = 0
  res@tmYLTickEndF        = 280; 250
  res@tmYLTickSpacingF    = 50;
  res@tmYLMinorPerMajor   = 4
  res@tmYLAutoPrecision   = True;False
  res@tmYLPrecision       = 3


  resR@tmYROn              = True
  resR@tmYRLabelsOn        = True
  resR@tiYAxisString = "Percent population (%)";+" [dashed]"
  resR@trYMinF         = 0
  resR@trYMaxF         := 40; 35.56
  resR@tmYRMode            = "Manual"
  resR@tmYRTickStartF      = 0
  resR@tmYRTickEndF        = 40; 36
  resR@tmYRTickSpacingF    = 5
  resR@tmYRMinorPerMajor   = 4
  resR@tmYRAutoPrecision   = False
  resR@tmYRPrecision       = 2

  res@xyLineColors     = "purple2";
  res@xyDashPatterns   = 0
  resR@xyLineColors    =  res@xyLineColors
  plot2(0) = gsn_csm_xy2(wks1, xx, G0_wpop_sun_median, G0_wpct_sun_median, res, resR);
  ;plot2(0) = gsn_csm_xy(wks1, xx, G0_wpct_sun_median, res); full solar
  gsres@gsFillColor       = "purple"                    ; color chosen
  dummy22 = gsn_add_polygon (wks1,plot2(0),xp,yp11,gsres);
  res@xyDashPatterns   = 1


  res@tiXAxisString    = " "
  res@tiYAxisString    = " "
  res@gsnLeftString    = " ";
  res@gsnCenterString  = " ";
  res@tmYLOn              = False
  res@tmYLLabelsOn        = False
  res@tmXBOn              = False
  res@tmXBLabelsOn        = False
  resR@tmYROn              = False
  resR@tmYRLabelsOn        = False
  resR@tiYAxisString      = ""

  res@xyLineColors     = "orange2";
  res@xyDashPatterns   = 0; 1
  resR@xyLineColors    =  res@xyLineColors
  plot2(1) = gsn_csm_xy2(wks1, xx, G0_wpop_diff_median, G0_wpct_diff_median, res, resR);
  ;plot2(1) = gsn_csm_xy(wks1, xx, G0_wpct_diff_median, res); diffuse
  gsres@gsFillColor       = "orange"                    ; color chosen
  dummy11 = gsn_add_polygon (wks1,plot2(1),xp,yp22,gsres);
  res@xyDashPatterns   = 1

  res@xyLineColors      = "slategray"
  res@xyDashPatterns    = 0; 
  resR@xyLineColors    =  res@xyLineColors
  plot2(2) = gsn_csm_xy2(wks1, xx, G0_wpop_zero_median, G0_wpct_zero_median, res, resR);
  ;plot2(2) = gsn_csm_xy(wks1, xx, G0_wpct_zero_median, res);
  gsres@gsFillColor       = "slategray4"               ; color chosen
  dummy33 = gsn_add_polygon (wks1,plot2(2),xp,yp33,gsres);


  ;===add error bars (+-50% of Rd for Outdoor shade scenario)
  ;gsn_add* templates are functions that we set to dummy values. Since
  ;we are going to draw numerous error bars, we create two arrays to hold the dummy values.
  error_bar1 = new(nlev,graphic)
  error_bar2 = new(nlev,graphic)
 
  xp2    := new( (/2,nlev/), float )
  xp2(0,:) = xx
  xp2(1,:) = xx
  ;diffuse scenario: add range of 10%-60% Rs  
  yp3_area  = new( (/2,nlev/), float )
  yp3_pop   = new( (/2,nlev/), float )
  yp3_area(0,:) = G0_warea_diff2_median ; uppper bound
  yp3_area(1,:) = G0_warea_diff1_median ; lower bound
  yp3_pop(0,:) = G0_wpop_diff2_median ; uppper bound
  yp3_pop(1,:) = G0_wpop_diff1_median ; lower bound

  res3                   = True                    ; res for confidence interval
  res3@gsLineThicknessF     = 1
  ;res3@xyCurveDrawOrder = "PostDraw"
  res3@trYAxisType = "LinearAxis"; linear scale for error bars
 
  do lev=0,nlev-1
    ;res3@gsLineDashPattern  = 2
    res3@gsLineColor        = "orange2";
    error_bar1(lev) = gsn_add_polyline(wks1,plot1(1),xp2(:,lev),yp3_area(:,lev), res3)
    error_bar2(lev) = gsn_add_polyline(wks1,plot2(1),xp2(:,lev),yp3_pop(:,lev), res3)
  end do

  ;--add some text
  txres                     = True                 ; text mods desired
  txres@txFontHeightF       = 0.012                ; default size is HUGE!
  txres@txJust              = "CenterLeft"         ; puts text on top of bars
  ;gsn_text(wks1,plot1(2),"Global scale", 1.2,900,txres) ;
  ;gsn_text(wks1,plot2(2),"*masked by rural population data", 1.2,res@trYMaxF*0.95,txres) ; on top


  ;--add legend for shade
  lbres                    = True
  lbres@lbBoxMajorExtentF  = 0.6           ; smaller value for more space between color boxes
  lbres@lbMonoFillPattern  = True
  lbres@lbPerimOn          = False
  lbres@lbBoxLinesOn       = False
  lbres@lbFillOpacityF     = 0.2; 0.1

  ;legend for lines, more general function gsn_legend_ndc
  lgres                    = True
  lgres@lgLineThicknessF   = 1.5                 ; line thicknesses
  lgres@lgPerimOn          = False               ; no perimeter
  lgres@lgLabelOffsetF     = 0.08
  lgres@lgLabelFontHeightF = 0.07;          ; font height (also scales with vpWidthF)

  ncase=3; 2
  if (ncase .eq. 3) then
  ;--legend for lines, more general function gsn_legend_ndc
  lbres@vpWidthF           = 0.18          ; labelbar width
  lbres@vpHeightF          = 0.08          ; labelbar height
  labels := (/"","",""/); no labels for the shade bars
  lbres@lbFillColors       := (/"slategray4","orange","purple"/) 
  gsn_labelbar_ndc(wks1,3,labels,0.07,0.9,lbres) ; shade

  lgres@vpWidthF           = 0.11                ; width of legend (NDC)
  lgres@vpHeightF          = 0.08           ; height of legend (NDC)
  lgres@lgDashIndexes      := (/0,0,0/)          ; line patterns
  lgres@lgLabelFontHeightF = 0.07;          ; font height (also scales with vpWidthF)
  lgres@lgLineColors       := (/"slategray","orange2","purple2"/)
  labels := (/"Dark","Shade","Sun"/)
  gsn_legend_ndc(wks1,3,labels,0.12,0.9,lgres)

  else

  ;--legend for lines, more general function gsn_legend_ndc
  lbres@vpWidthF           = 0.18          ; labelbar width
  lbres@vpHeightF          = 0.06;0.08          ; labelbar height
  labels := (/"",""/); no labels for the shade bars
  lbres@lbFillColors       := (/"slategray4","purple"/)
  gsn_labelbar_ndc(wks1,2,labels,0.07,0.9,lbres) ; shade

  lgres@vpWidthF           = 0.1; 0.11                ; width of legend (NDC)
  lgres@vpHeightF          = 0.06;0.08           ; height of legend (NDC)
  lgres@lgDashIndexes      = (/0,0,0/)          ; line patterns
  lgres@lgLabelFontHeightF = 0.07;          ; font height (also scales with vpWidthF)
  lgres@lgLineColors       := (/"slategray","purple2"/)
  labels := (/"Dark","Sun"/)
  gsn_legend_ndc(wks1,2,labels,0.12,0.9,lgres)

  end if

  print("finish plot b")

  draw(plot1)
  draw(plot2)



  ;--------Fig4cd areal and population impacts (agricultural workers)------------

  plot3 = new(3,graphic)
  plot4 = new(3,graphic)

  res@vpYF          = 0.50   ; location of plot1
  res@vpXF          = 0.10    ; location of plot1

  res@tmYLOn              = True
  res@tmYLLabelsOn        = True
  res@tmXBOn              = True
  res@tmXBLabelsOn        = True
  res@trYAxisType = "LinearAxis"
  res@gsnStringFontHeightF = 0.014;
  res@gsnLeftString       = "~F22~c";
  res@gsnCenterString     = "Focus on agricultural workers"; "Global";
  res@gsnLeftStringOrthogonalPosF = 0.03 ;adjust the space/distance of gsnLeftString/gsnCenterString
  res@gsnCenterStringOrthogonalPosF = 0.02 ;adjust the space/distance of gsnLeftString/gsnCenterString
  ;res@tiYAxisString    = "Million hectares (Mha)"; "  Land area (Mha)~C~       with G~S~day~N~>0"
  ;res@tiYAxisString    = "Land area with G~S~day~N~>0 or T~B~W~N~~S~day~N~>35 ~C~                       (Mha)"; 
  res@tiYAxisString    = "Land area with G~S~day~N~>0 (Mha)"; 
  res@tiYAxisJust      = "CenterCenter"
  res@tiYAxisOffsetXF  = -0.01; -0.01
  res@tiXAxisString    = "Warming since pre-industrial (~S~o~N~C)"
  res@trYMinF             = 0; divide plot to a half
  res@trYMaxF             = 2628; 2190; Mha
  res@tmYROn              = False
  res@tmYLMode            = "Manual"
  res@tmYLTickStartF      = 0
  res@tmYLTickEndF        = 2600; 2100
  res@tmYLTickSpacingF    = 500
  res@tmYLMinorPerMajor   = 4
  res@tmYLAutoPrecision   = False
  res@tmYLPrecision       = 4

 ;resources for "right" Y axis
  resR = True
  resR@vpYF          = 0.50   ; location of plot1
  resR@vpXF          = 0.10    ; location of plot1
  resR@tiYAxisString = "Percent area (%)";+" [dashed]"
  resR@trYMinF         = 0
  resR@trYMaxF         = 30; 25
  resR@xyDashPatterns    = 0; 
  resR@xyLineThicknesses  = 1.5
  resR@xyLineOpacities    = 1
  resR@pmLegendDisplayMode    = "Never"

  resR@tmYRMode            = "Manual"
  resR@tmYRTickStartF      = 0
  resR@tmYRTickEndF        = 30; 25
  resR@tmYRTickSpacingF    = 5
  resR@tmYRMinorPerMajor   = 4
  resR@tmYRAutoPrecision   = False
  resR@tmYRPrecision       = 2
  resR@tmYROn              = True
  resR@tmYRLabelsOn        = True

  res@xyLineColors     = "purple2";
  resR@xyLineColors    =  res@xyLineColors
  res@xyDashPatterns   = 0
  ;plot3(0) = gsn_csm_xy(wks1, xx, G0_rarea_sun_median, res); full solar
  plot3(0) = gsn_csm_xy2(wks1, xx, G0_rarea_sun_median, G0_rareafrc_sun_median, res, resR); 
 

  gsres                   = True                        ; poly res
  gsres@gsFillOpacityF    = 0.1
  gsres@gsFillColor       = "purple"                    ; color chosen
  dummy4 = gsn_add_polygon (wks1,plot3(0),xp,yp4,gsres);
  res@xyDashPatterns   = 1


  res@tiXAxisString    = " "
  res@tiYAxisString    = " "
  res@gsnLeftString    = " ";
  res@gsnCenterString    = " ";
  res@tmXBOn              = False
  res@tmYLOn              = False
  res@tmXBLabelsOn        = False
  res@tmYLLabelsOn        = False
  resR@tmYROn             = False
  resR@tmYRLabelsOn       = False
  resR@tiYAxisString      = ""

  res@xyLineColors     = "orange2";
  res@xyDashPatterns   = 0; 1
  resR@xyLineColors    =  res@xyLineColors

  ;plot3(1) = gsn_csm_xy(wks1, xx, G0_rarea_diff_median, res); diffuse
  plot3(1) = gsn_csm_xy2(wks1, xx, G0_rarea_diff_median, G0_rareafrc_diff_median, res, resR); 
  gsres@gsFillColor       = "orange"                    ; color chosen
  dummy5 = gsn_add_polygon (wks1,plot3(1),xp,yp5,gsres);
  res@xyDashPatterns   = 1

  res@xyLineColors      = "slategray"
  resR@xyLineColors    =  res@xyLineColors
  res@xyDashPatterns    = 0; 3
  ;plot3(2) = gsn_csm_xy(wks1, xx, G0_rarea_zero_median, res);
  plot3(2) = gsn_csm_xy2(wks1, xx, G0_rarea_zero_median, G0_rareafrc_zero_median, res, resR); 
  gsres@gsFillColor       = "slategray4"               ; color chosen
  dummy6 = gsn_add_polygon (wks1,plot3(2),xp,yp6,gsres);

  print("finish plot c")



 ;use vpXF to place two plots horizontally
  res@vpXF          = 0.6       ; location of plot4
  resR@vpXF          = 0.6       ; location of plot4
  res@tmYLOn              = True
  res@tmYLLabelsOn        = True
  res@tmXBOn              = True
  res@tmXBLabelsOn        = True

  res@gsnLeftString       = "~F22~d";
  res@gsnCenterString     = "Focus on agricultural workers"; "Global";
  res@tiXAxisString    = "Warming since pre-industrial (~S~o~N~C)"
  ;res@tiYAxisString    = "Population with G~S~day~N~>0 or T~B~W~N~~S~day~N~>35 ~C~                      (million)"; 
  res@tiYAxisString    = "Population with G~S~day~N~>0 (million)"; 
  res@tiYAxisJust      = "CenterCenter"
  res@tiYAxisOffsetXF  = 0;0.01
  res@trYMinF             = 0; 
  res@trYMaxF             = 255; 212; million
  res@tmYLOn              = True
  res@tmYROn              = False
  res@tmYLLabelsOn        = True
  res@tmYRLabelsOn        = False
  res@tmYLMode            = "Manual"
  res@tmYLTickStartF      = 0
  res@tmYLTickEndF        = 250; 210
  res@tmYLTickSpacingF    = 50;
  res@tmYLMinorPerMajor   = 4
  res@tmYLAutoPrecision   = True;False
  res@tmYLPrecision       = 3


  resR@tiYAxisString = "Percent population (%)";+" [dashed]"
  resR@tmYRMode            = "Manual"
  resR@tmYRTickStartF      = 0
  resR@tmYRTickEndF        = 30; 25
  resR@tmYRTickSpacingF    = 5
  resR@tmYRMinorPerMajor   = 4
  resR@tmYRAutoPrecision   = False
  resR@tmYRPrecision       = 2
  resR@tmYROn              = True
  resR@tmYRLabelsOn        = True

  res@xyLineColors     = "purple2";
  res@xyDashPatterns   = 0
  resR@xyLineColors    =  res@xyLineColors
  plot4(0) = gsn_csm_xy2(wks1, xx, G0_rpop_sun_median, G0_rpct_sun_median, res, resR);
  ;plot4(0) = gsn_csm_xy(wks1, xx, G0_rpct_sun_median, res); full solar
  gsres@gsFillColor       = "purple"                    ; color chosen
  dummy44 = gsn_add_polygon (wks1,plot4(0),xp,yp44,gsres);
  res@xyDashPatterns   = 1


  res@tiXAxisString    = " "
  res@tiYAxisString    = " "
  res@gsnLeftString    = " ";
  res@gsnCenterString  = " ";
  res@tmXBOn              = False
  res@tmYLOn              = False
  res@tmXBLabelsOn        = False
  res@tmYLLabelsOn        = False
  resR@tmYROn             = False
  resR@tmYRLabelsOn       = False
  resR@tiYAxisString      = ""

  res@xyLineColors     = "orange2";
  res@xyDashPatterns   = 0; 1
  resR@xyLineColors    =  res@xyLineColors
  plot4(1) = gsn_csm_xy2(wks1, xx, G0_rpop_diff_median, G0_rpct_diff_median, res, resR);
  ;plot4(1) = gsn_csm_xy(wks1, xx, G0_rpct_diff_median, res); diffuse
  gsres@gsFillColor       = "orange"                    ; color chosen
  dummy55 = gsn_add_polygon (wks1,plot4(1),xp,yp55,gsres);
  res@xyDashPatterns   = 1

  res@xyLineColors      = "slategray"
  res@xyDashPatterns    = 0; 
  resR@xyLineColors    =  res@xyLineColors
  plot4(2) = gsn_csm_xy2(wks1, xx, G0_rpop_zero_median, G0_rpct_zero_median, res, resR);
  ;plot4(2) = gsn_csm_xy(wks1, xx, G0_rpct_zero_median, res);
  gsres@gsFillColor       = "slategray4"               ; color chosen
  dummy66 = gsn_add_polygon (wks1,plot4(2),xp,yp66,gsres);


  ;===add error bars (+-50% of Rd for Outdoor shade scenario)
  ;gsn_add* templates are functions that we set to dummy values. Since
  ;we are going to draw numerous error bars, we create two arrays to hold the dummy values.
  error_bar1 = new(nlev,graphic)
  error_bar2 = new(nlev,graphic)
 
  xp2    := new( (/2,nlev/), float )
  xp2(0,:) = xx
  xp2(1,:) = xx
  ;diffuse scenario: add range of 10%-60% Rs  
  yp3_area  = new( (/2,nlev/), float )
  yp3_pop   = new( (/2,nlev/), float )
  yp3_area(0,:) = G0_rarea_diff2_median ; uppper bound
  yp3_area(1,:) = G0_rarea_diff1_median ; lower bound
  yp3_pop(0,:) = G0_rpop_diff2_median ; uppper bound
  yp3_pop(1,:) = G0_rpop_diff1_median ; lower bound

  res3                   = True                    ; res for confidence interval
  res3@gsLineThicknessF     = 1
  ;res3@xyCurveDrawOrder = "PostDraw"
  res3@trYAxisType = "LinearAxis"; linear scale for error bars
 
  do lev=0,nlev-1
    ;res3@gsLineDashPattern  = 2
    res3@gsLineColor        = "orange2";
    error_bar1(lev) = gsn_add_polyline(wks1,plot3(1),xp2(:,lev),yp3_area(:,lev), res3)
    error_bar2(lev) = gsn_add_polyline(wks1,plot4(1),xp2(:,lev),yp3_pop(:,lev), res3)
  end do


  ;--add some text
  txres                     = True                 ; text mods desired
  txres@txFontHeightF       = 0.012                ; default size is HUGE!
  txres@txJust              = "CenterLeft"         ; puts text on top of bars
  ;gsn_text(wks1,plot3(2),"Global scale", 1.2,900,txres) ;
  ;gsn_text(wks1,plot4(2),"*masked by rural population data", 1.2,res@trYMaxF*0.95,txres) ; on top


  ;--add legend for shade
  lbres                    = True
  lbres@lbBoxMajorExtentF  = 0.6           ; smaller value for more space between color boxes
  lbres@lbMonoFillPattern  = True
  lbres@lbPerimOn          = False
  lbres@lbBoxLinesOn       = False
  lbres@lbFillOpacityF     = 0.2; 0.1

  ;legend for lines, more general function gsn_legend_ndc
  lgres                    = True
  lgres@lgLineThicknessF   = 1.5                 ; line thicknesses
  lgres@lgPerimOn          = False               ; no perimeter
  lgres@lgLabelOffsetF     = 0.08
  lgres@lgLabelFontHeightF = 0.07;          ; font height (also scales with vpWidthF)

;  ncase=3; 2
;  if (ncase .eq. 3) then
;  ;--legend for lines, more general function gsn_legend_ndc
;  lbres@vpWidthF           = 0.18          ; labelbar width
;  lbres@vpHeightF          = 0.08          ; labelbar height
;  labels := (/"","",""/); no labels for the shade bars
;  lbres@lbFillColors       := (/"slategray4","orange","purple"/) 
;  gsn_labelbar_ndc(wks1,3,labels,0.07,0.5,lbres) ; shade
;
;  lgres@vpWidthF           = 0.11                ; width of legend (NDC)
;  lgres@vpHeightF          = 0.08           ; height of legend (NDC)
;  lgres@lgDashIndexes      := (/0,0,0/)          ; line patterns
;  lgres@lgLabelFontHeightF = 0.07;          ; font height (also scales with vpWidthF)
;  lgres@lgLineColors       := (/"slategray","orange2","purple2"/)
;  labels := (/"Dark","Shade","Sun"/)
;  gsn_legend_ndc(wks1,3,labels,0.12,0.5,lgres)
;
;  else
;
;  ;--legend for lines, more general function gsn_legend_ndc
;  lbres@vpWidthF           = 0.18          ; labelbar width
;  lbres@vpHeightF          = 0.06;0.08          ; labelbar height
;  labels := (/"",""/); no labels for the shade bars
;  lbres@lbFillColors       := (/"slategray4","purple"/)
;  gsn_labelbar_ndc(wks1,2,labels,0.07,0.5,lbres) ; shade
;
;  lgres@vpWidthF           = 0.1; 0.11                ; width of legend (NDC)
;  lgres@vpHeightF          = 0.06;0.08           ; height of legend (NDC)
;  lgres@lgDashIndexes      = (/0,0,0/)          ; line patterns
;  lgres@lgLabelFontHeightF = 0.07;          ; font height (also scales with vpWidthF)
;  lgres@lgLineColors       := (/"slategray","purple2"/)
;  labels := (/"Dark","Sun"/)
;  gsn_legend_ndc(wks1,2,labels,0.12,0.5,lgres)
;
;  end if

  print("finish plot d")


  draw(plot3)
  draw(plot4)
  frame(wks1)




end 
